################################################################
##### Calculating WTP from SEP Wave 8 Conjoint Data  ######
################################################################

rm(list = ls())


# Load required packages
library(cregg)
library(ggplot2)
library(dfidx)
library(mlogit)
library(dplyr)
library(readr)

# Load conjoint data
dfconj <- read_csv("/Users/cbrugge/Desktop/CE Paper/wtp_data_ce.csv") %>% 
  as.data.frame() %>% 
  select(-"...1")


# Prep data for dfidx  ----------------------------------------------------

dfconj_prod1 <- subset(dfconj, product == 1)

colnames(dfconj_prod1) <- c("ID", "product" , "roundid", "choice", "cost.A", "life.A", "env.A", 
                           "eff.A", "recy.A", "rep.A", "product.type")

dfconj_prod1 <- subset(dfconj_prod1, select = -c(product))

dfconj_prod2 <- subset(dfconj, product == 2)

colnames(dfconj_prod2) <- c("ID", "product" , "roundid", "choice", "cost.B", "life.B", "env.B", 
                           "eff.B", "recy.B", "rep.B", "product.type")


dfconj_prod2 <- subset(dfconj_prod2, select = -c(product, choice, roundid, ID, product.type))

data <- cbind(dfconj_prod1, dfconj_prod2)

col_order <- c("ID", "roundid", "choice", "cost.A", "life.A", "env.A", "eff.A", "recy.A", "rep.A", 
               "cost.B", "life.B", "env.B", "eff.B", "recy.B", "rep.B", "product.type")

data <- data[, col_order]

# convert relevant variables to characters
data <- data %>%
  mutate(ID = as.character(ID),
         roundid = as.character(roundid),
         choice = as.character(choice))

#add character string to ID 
data$ID <- paste("R_", data$ID, sep = "")

data <- data %>% mutate(choice = ifelse(choice == "1", "A", "B"))

# Create a new unique identifier by combining "round_id" and "ID"
data$unique_id <- as.character(interaction(data$roundid, data$ID, sep = "."))


# create individual datasets by product -----------------------------------
data_smartphone <- subset(data, product.type == "Smartphone")
data_tv<- subset(data, product.type == "TV")
data_washingmachine<- subset(data, product.type == "Washing machine")

# # Convert them to the new values 250, 850, 1900, and 3100 to allow non-linear analysis
data_tv <- data_tv %>%
  mutate(cost.A = case_when(
    cost.A == 1 ~ 250,
    cost.A == 2 ~ 850,
    cost.A == 3 ~ 1900,
    cost.A == 4 ~ 3100,
    TRUE ~ cost.A  
  ),
  cost.B = case_when(
    cost.B == 1 ~ 250,
    cost.B == 2 ~ 850,
    cost.B == 3 ~ 1900,
    cost.B == 4 ~ 3100,
    TRUE ~ cost.B  
  ))

# apply dfidx --------------------------------------------------------------

#mlogdata <- dfidx(data, idx = "unique_id", choice = "choice", varying = 4:15, sep = ".")
mlogdata <- dfidx(data_tv, idx = "unique_id", choice = "choice", varying = 4:15, sep = ".")

# calculate mlogit and wtp ------------------------------------------------
m1_tv <- mlogit(choice ~ cost + life + env + eff + recy + rep | -1,
             mlogdata)
summary(m1_tv)

# Print a formatted regression table
stargazer(m1_tv, type = "text", covariate.labels = c("Price","Function duration","Environmental impact during production",
                                                  "Energy efficiency","Recyclability","Repairability"), 
          dep.var.labels = "Product choice")
# Save the regression table as an HTML file
stargazer(m1_tv, type = "html", covariate.labels = c("Price","Function duration","Environmental impact during production",
                                                  "Energy efficiency","Recyclability","Repairability"), dep.var.labels =
            "Product choice (TV)", out = "/Users/cbrugge/Desktop/CE Paper/mixed_logit_tv.html")

# odds ratios
odds <- exp(coef(m1_tv, matrix = TRUE))
precentage_increase_by_level <- 1-odds
precentage_increase_by_level


wtp <- coef(m1_tv)[2:length(coef(m1_tv))] / coef(m1_tv)["cost"]
wtp
#After running this code, the wtp variable will contain the calculated WTP values
#for each attribute level (excluding the intercept) relative to a one-point increase
#in the price attribute.

# Transform the WTP values to percentage change
wtp_percentage_change <- (exp(wtp) - 1) * 100
print(wtp_percentage_change)

